Soliton-radiation coupling in the parametrically driven, damped 

nonlinear Schrodinger equation 

V. S. ShchesnovicbQ and I. V. BarashenkovQ 

Department of Mathematics and Applied Mathematics, University of 
Cape Town, Private Bag Rondebosch 7701, South Africa. 



1 valery @ mat hs.uct.ac.za 
2 igor@cenerentola.mth.uct.ac.za 



Abstract 



We use the Riemann-Hilbert problem to study the interaction of the soliton with 
radiation in the parametrically driven, damped nonlinear Schrodinger equation. The 
analysis is reduced to the study of a finite-dimensional dynamical system for the am- 
plitude and phase of the soliton and the complex amplitude of the long-wavelength 
radiation. In contrast to previously utilised Inverse Scattering-based perturbation 
techniques, our approach is valid for arbitrarily large driving strengths and damping 
coefficients. We show that, contrary to suggestions made in literature, the complexity 
observed in the soliton's dynamics cannot be accounted for just by its coupling to the 
long-wavelength radiation. 



1 Introduction 



A variety of nonlinear wave phenomena in one dimension can be modelled by the perturbed 
nonlinear Schrodinger (NLS) equation: 



iqt + q xx + 2\q\ 2 q = K{t, x, q, q*, q x , q* x , ...). 



(1) 



(Here and below the asterisk denotes complex conjugation.) The present paper deals with 
the parametrically driven, damped NLS, for which 



TZ = he 2inot q* - i^q; h, Q Q , 7 > 0. 



(2) 



This equation describes nonlinear Faraday resonance in a vertically oscillating water trough [I], 
0, |3], £|, |BJ, H, |7|, |8|, ||, [L(], [11], |12| ; an easy-plane ferromagnet with a combination of a stationary 



and a high-frequency magnetic field in the easy plane fll3[ ; and the effect of phase-sensitive 
amplifiers on solitons propagating in optical fibres fl4| , [15 , [16 . 

The equation ([l])-(|2]) has two stationary soliton solutions, one of which is unstable for 
all h and 7 and hence usually disregarded fI5||. (The frequency fl can always be scaled 



to unity leaving h and 7 as the only two control parameters.) The other stationary soliton 
is stable for small values of h but looses its stability to a periodically oscillating soliton as 
h is increased above a certain critical value h = h c (j), for the fixed 7 The chart of 

attractors arising as h is increased further, was compiled in [17| and can be summarised 
as follows. If 7 is small, the oscillating soliton undergoes a sequence of bifurcations which 
culminate in a soliton whose amplitude, width and phase are changing chaotically in time. 
For an even greater h (and the same, fixed, 7) the soliton breaks up and decays to zero 
whereas increasing the driver's strength still further, a spatio-temporal chaotic state sets in. 
On the other hand, if 7 is large, the bifurcations of the soliton's period and its breakup are 
not observed; instead, as h is increased for this 7, the periodically oscillating soliton yields 
directly to the spatio-temporal chaos. 

Some insight into the mechanism controlling the soliton's transformations was gained in 
Ref . [T8[] where a reduced amplitude equation was derived for the perturbation of the station- 
ary soliton. Its analysis demonstrated that the emission of radiation waves plays a major 
role in the soliton's dynamics. In particular, the soliton-radiation interaction accounts for 
the stable periodic oscillations of the damped soliton (7 7^ 0) and for the absence of the 
periodicity in the undamped situation, 7 = 0. However, the analysis of [|l8j was confined 
to a neighbourhood of the instability threshold h = h^) and hence the reduced ampli- 
tude equation could not reproduce the entire complexity in the soliton's dynamics which is 
observed in numerical simulations with larger h |T7] . 

In the present work we study the soliton-radiation interaction from a different perspec- 
tive. Our decomposition of the phase space into soliton and radiation modes will be based 
not on the linearisation about the stationary soliton in the "g-space" (which was the ap- 
proach of Ref. [pl|), but on the analysis of the spectral data of the Riemann-Hilbert problem 
associated with the unperturbed NLS equation. This will allow us to examine the role of 
the soliton-radiation coupling for arbitrary h and 7, and not just in the neighbourhood of 
the instability onset. 



There is a number of perturbation schemes available in literature which exploit the 
proximity of the perturbed NLS (|]) to the completely integrable case of the "pure" NLS, 
eq.(0) with TZ = 0. (See g§ fZ5[ 0, [2§, 0, and a review |^). Solutions 
of equation (JJ) have their images in the space of spectral data for any TZ, of course, but 
this is of little use in the general case. The difficulty here is that the evolution equations 
for the spectral data involve the associated eigenfunctions (or, equivalently, solutions of 
the Riemann-Hilbert problem). In order to obtain a closed system of equations for the 
spectral data, one assumes that TZ is small in some sense and expands the spectral data and 
eigenfunctions in powers of the small parameter [p4j |25fl . In the adiabatic approximation, 
for example, one first derives equations of the (slow) evolution of the soliton's parameters 
(ignoring the radiation completely) and then calculates the spectral density of the radiation 
emitted by this soliton |2^, The back-reaction of the radiation on the soliton is 
not taken into account in this approach, therefore. The adiabatic approximation is capable 
of capturing some basic essentials of the soliton's dynamics, such as the phase-locking of 
the soliton to the periodic driver (see e.g. [^6|) or stability against perturbations of its 
parameters f25|, and is usually sufficient for very small right-hand sides in ([[]). (See |32| for 



details.) However, it becomes inadequate for somewhat larger TZ, where the back-reaction 
of the emitted radiation on the soliton cannot be disregarded [f28| , |29| , j33fl . (For example, the 
adiabatic approximation does not capture the oscillatory instability of the parametrically 
driven soliton which sets in for 7 = and h as small as 0.064 |l3f .) 

An attempt to go beyond the adiabatic approximation and consider the radiation de- 
grees of freedom on equal footing with those of the soliton was made by the authors of 
Ref. |33j whose finite-dimensional reduction included the complex amplitude of the k = 0- 
radiation. However, their derivation was not entirely self-consistent; in particular, their 
approach produced an equation for the radiation part of the spectral data which did not 
contain a damping term. To reconcile conclusions of their finite-dimensional analysis with 
direct numerical simulations of the full partial differential equation, the authors had to add 
the damping in an ad hoc way. The amplitude of the radiation wave also remained undefined 



and had to be chosen so as to match the numerics 33 



Similarly to Ref. |33fl, the purpose of the present work is to study the effect of the 
soliton-radiation coupling on the internal dynamics of the damped-driven soliton. However, 
unlike the analysis of Ref. []33 and perturbation schemes appeared elsewhere, our approach 
is not using the assumption of the smallness of TZ. Instead, we will exploit the fact that 
the stationary soliton of equation (|l])-(|]) with TZ ^ coincides — up to a simple phase 
transformation — with the soliton solution of equation ([!]) with TZ = 0. Consequently, 
the stationary soliton of the parametrically driven damped NLS with arbitrarily large h 
and 7 can be associated with a stationary zero of the Riemann-Hilbert problem underlying 
the "pure", integrable, NLS equation. This observation allows to choose a different small 
parameter; instead of the smallness of h and 7 we will be utilising the proximity of the 
solution in question to the stationary soliton of the perturbed equation. (Another assumption 
that we are going to make, following is that the radiation is linear, i.e. it couples to the 
nonlinearly evolving soliton but does not interact with itself.) Using this small parameter we 
will be able to obtain a closed system of equations describing, approximately, the evolution 
of the spectral data. 



Some analytical insights can be gained already from the linearisation of this system 
(about the zero of the Riemann-Hilbert problem corresponding to a single soliton). We will 
show that the linearisation explains the origin of the oscillatory instability which serves as 
a starting point of the sequence of secondary bifurcations and leads to the emergence of 
the increasing complexity in the soliton's dynamics. (So far, the oscillatory instability and 
related Hopf bifurcation remained just facts of numerical analysis Ul3|, ||18|| .) It also indicates 
that the soliton interacts most intensively with radiation waves near the lower boundary 
of their spectrum (where k = 0). This seems to be in agreement with earlier suggestions 
- made for a closely related externally driven NLS — that keeping just the k = mode 
is sufficient to capture the basic features of the partial differential equation |33|, [34], |35|j . 
Accordingly, we focus our subsequent efforts on the verification of this hypothesis. Namely, 
we explore the effect of the coupling to the k = radiation on the nonlinear dynamics of 
the soliton. Keeping only infinitely long waves allows to obtain a four- dimensional system 
for the soliton phase and the complex amplitude of the radiation. Results of the analysis of 
this system are then compared to the phenomenology of the parametrically driven soliton 
reported in literature. We will show that taking the k = radiation into account can explain 
some dynamical effects, most notably the occurrence of the soliton's breakup and decay to 
zero. We will also identify aspects of the behaviour that cannot be attributed just to the 
coupling of the soliton to the long- wavelength radiation and therefore require invoking other 
degrees of freedom. These aspects will include, in particular, the shape of the instability 
domain on the (h, 7)-plane and the route to the (temporal) chaos. 

The paper is organised as follows. The next section contains a brief summary of the 
Riemann-Hilbert problem and the perturbation theory based upon it. Section || discusses 
the adiabatic approximation and its shortcomings. In section [| we linearise the evolution 
equations for the spectral data, while the nonlinear four-dimensional system for the soliton's 
parameters and the complex amplitude of the k = radiation is derived in section [5]. 
Numerical simulations of this system are reported in section |6|. The last section (section [7]) 
contains conclusions of our analysis. 



2 Inverse Spectral Transfrom for the "pure" and per- 
turbed NLS equation 

In this section we review the main points of the modern version of the Inverse Scattering 
Transform, known as the method of the Riemann-Hilbert problem, for the NLS equation. 
After that we outline the basic principles of the Inverse Scattering-based perturbation theory, 
in its particular Riemann-Hilbert formulation. 



2.1 The Riemann-Hilbert problem 

The applicability of the Inverse Scattering Transform to the unperturbed NLS equation (eq. 
([[]) with 1Z — 0) is due to the fact that the NLS serves as the compatibility condition for 
the following system of two linear equations for the matrix-valued function ty(x, t; (): 

d x ^ = iC[a 3 ,^]+iQ^, (3) 



dtV = -2i( 2 [a 3 , *] - (2<g + a 3 Q x - iQ 2 a z ) (4) 
where the potential matrix 

Q=( ° q 

g = q(x,t) is a solution of the unperturbed NLS and er 3 is the Pauli matrix |Tj|. Knowing 
the potential Q can be recovered from the asymptotic expansion of ^(C) as £ ^ oo: 

*(C)=/ + vI>«C- 1 + ..., Q=[*M,a a ], (5) 

whereas ^ is found via the analytic factorisation in the complex £-plane. The factorisation 
problem is known as the matrix Riemann-Hilbert problem; the interested reader may consult 
Refs. |2T], ^2], |2Bj for the full account of this technique while here we only give a brief 
summary. 

First one defines the Jost solutions of equation ([$]) by their asymptotic behaviour as 
x — > ±oo: J± — > I. Then the matrix-valued function 

*+(C) = ( (J + HO {J-UO) (6) 

is a solution to (y) - ® holomorphic in the upper half of the complex £-plane (Im£ > 0). In 
®j (J±)-i stands for the Z-th column of J±. Noting that the linear problem admits 
an involution ^(C) — ^ X (C*)? we introduce a matrix function holomorphic in the lower 
half-plane: 

*= 1 (C) = *l(O = ((^: 1 )i.(0, U- 1 ) 2 -(C)) T , (7) 

where ( J±~ 1 )i- denotes the /-th row of the matrix J±~ 1 and T indicates transposition. The 
functions ^+(C) an d ^-(C), solutions to equation @, can be expressed through the Jost 
solutions and the elements of the scattering matrix, defined as 

S = e" iCxCT3 J" 1 J_e iCxa3 . (8) 

Introducing the upper and lower-triangular matrices S±, satisfying the equation S + = SS-, 
by 

S+= {1 £)' S - = {[s- 1 hi i)' (9) 

we have 

^ + = ./.,'-'""•. S\ f ^I 1 = e^^tg-iC^a j-i ( 10 ) 

This leads to the Riemann-Hilbert problem of finding the matrix- valued functions X I / +(C) 
and ^/I 1 (C), holomorphic in the upper and lower half-plane of (, respectively, and satisfying 

^Z 1 ^+ = e iixas Ge' iixa3 , (11) 

on the real line |21||. Here *f>± — > I as £ — ► oo and 



G 



1 9 

9* 1 



with g(C) = Si2(C)- The t-dependence of g follows from equations ( |TTD and (Q). We have 
G t = -2i( 2 [cr 3 , G], hence g t = -4i( 2 g. 

If the det \&+(C) has zeros in its analyticity domain, the Riemann-Hilbert problem is said 
to be singular (or with zeros). Owing to the involution, zeros of det\l/I 1 (C) are complex 
conjugates of those of det The solution to the Riemann-Hilbert problem with zeros 

can be written as 

*±(C) = *o±(C)r(C), (12) 

where det^o±(C) 7^ and the rational matrix function T(C) (the "dressing matrix") has 
zeros of det\l/ + (C) and poles at the zeros of det StC 1 ^). In the case of N simple zeros Q, 
j = 1, . . . , N (i.e. det \P+(C) — C(C — Cj) as C ~ > (j), the dressing matrix has the following 
structure (see also pi]): 

no^-f 'y^i as) 

where D ln = (( n - Cf) -1 (pj I Pn) and (pi \ p n ) = (p*)i(p n )i + (Pih&nh- The vector-columns 
| pj) and vector-rows (pj | are defined by 

9 + (C i )\p J )=0, (Pj | Vz\(j) = 0. (14) 

The involution property gives (pj \ = \ pj)^ and r _1 (C) = rt(£*). 

The matrix functions \&o±(C) (0) solve the following regular Riemann-Hilbert problem: 

^ol 1 ^ = iv— 67 0v! -T ; . ImC = 0, (15) 

where \I/o± — > I as C, — > oo. The solution to the regular Riemann-Hilbert problem is unique 
due to the normalisation condition at infinity. 

The spectral data defining a unique solution to the Riemann-Hilbert problem consists of 
two parts: the discrete set of Q and \p 3 ) (j = 1, N) and the continuous data g(Q = G 12(C). 
The pure soliton solutions arise from the Riemann-Hilbert problem with zeros provided 
g(C) = 0, i.e., * 0± (C) = I- 

To derive the coordinate dependence of the discrete data we note that det\l/±(C) is 
t- independent, hence (Q)t = 0. The coordinate dependence of | pj) is obtained by the 
differentiation of the first relation in flUD and using equations (|3|), (|) and (|T^): 

I Pj)x = I Pj), I Pj)t = -2< ? 2 (T3 I pj). (16) 



(0) 1 exp(a i + z% + 
Pj) = exp(/ i a 3 ) I ; ) = [ e _ fj ) , (17) 



Hence 

-/: 

where /j = — 2i£jt and we have defined exp(aj + i9j) = (pf )x/ (pf )2, with aj and 0j 
real constants. 

In what follows we will need the one-soliton dressing matrix. For N = 1 equations (13) 
and (0) yield 

sechz / (C - Ci)e z + (C - Ci*)e- 2 -2^ 
lU 2(C-C?)V -2^e"^ (C - CiV + (C - Ci)e" z 1 ' ' 



where we have defined 



z = a + 2Re/i = —2r\{x — xq), 



if = 9 + 2Im/i = —-z + <£>o; 

rj 



(19) 



decomposed Ci = £ + ^7, and denoted a\ = a and #1 = 9. In ( |H?D we have introduced the 
position x and the core phase (p of the soliton, via 



x=o 2r] 



(a + 877ft) , y?o = ip 



9 + 2£x + 4(r7 2 - i 2 )t. 



X = XQ 



(20) 



Finally, the one-soliton solution to the (unperturbed) NLS equation parametrised by the 
Riemann-Hilbert data reads 



q s (x,t) = -2 lim ^12 (C) = 2iije* p aedxz. 



(21) 



Here 2rj, 4f and a give, respectively, the amplitude, velocity and initial position of the 
soliton. 9 is the initial phase at the point x = 0: 9 = <p\ x=t=0 . 



2.2 Evolution of the spectral data in the perturbed NLS equation 

When TZ 7^ in equation ([I]), the evolution of the associated spectral data becomes nonlinear 
and complicated. Here our analysis follows the lines of Refs. |23], [25], [27], [Z8[ [2J| [TI[ . 
To distinguish between the integrable and perturbation-induced t-dependence, we use the 
"variational derivative" notation. For instance, the perturbation (|J) can be written as 

U = he 2iaot q*-i<yq = & 

ot 

Introducing an off-diagonal matrix 

-TV 

and the linear functional of the perturbation 

X 

±00 



R= t n ) ( 22 ) 



(23) 



we can write the variational derivative of in the following form (see f3lj for details): 



St 

where 



6* + {x, C) = Q e Kxa m ^^ Q e -K^^ (24) 



n / . ^ _ / -T n (x,oo;C) T 12 (-oo, ar; C) \ 
i H a;,y ^ _r 21 (x,oc;C) T 22 (-oo,x;C) ) 



(In (^), (pi ) and below we are omitting the explicit t-dependence for notational conve- 



nience.) The evolution functional H(x, Q is meromorphic in the upper half-plane of ( and 
has simple poles at zeros of det \I/+(C) (which are assumed to be simple). The r.h.s. of (^|) 
describes the perturbation- induced evolution of and should be added to the r.h.s. of (||). 
In view of the involution ^(C) — ^ -1 (C*)> the evolution of ^Z 1 is given by the Hermitian 
conjugate of equation (|24|): 



S ^-[( x > Q. = e iCxa m\x, ne'^VzUx, C). (25) 
ot 

Relations (p^)-(^5|) yield equations for the perturbation-induced evolution of the com- 
plete set of the spectral data {g(()] Cj, a j, = !)•••? ^)} : 

^ = -resT 22 (0), (26) 

^( % + iBi) = T 2 r 2 cs) (0) - exp Li j dt(j - a, - ie\ T^(Q, (27) 

-4iCM0 + T 12 (C)+T 22 (C)^(C), (28) 



MO 



dt 

where T(£) = T(— oo, oo; () and 

T^(Q = lim {T(0 - (C - Ci) _1 resT(0)} 

is the regular part of T(£) at the pole £ = Q. In the derivation of (^)- (|28|) we used equations 
TJ), (|TJ), flT7|) and the identities 



resn(x,0) | Pj) = | Pj), ^ + (x,Q)e l ^resU(xX 1 ) = 0, (29) 

where resII(C J ) denotes the residue of IT(£) at Q and | pj) = exp(—i(jxa 3 ) \ pj(x)). 

Equations (f26|)-(|2"8|) are highly nonlinear since the right-hand sides are dependent on the 
function \I/+(C) which is itself to be constructed from the spectral data. However, these 
equations can be simplified by expanding in powers of a suitably chosen small parameter. 
In section 4 we will make a choice of the small parameter that will allow us to develop a 
rigorous approach to the soliton-radiation interaction. For completeness of the presentation, 
however, we consider the standard adiabatic approximation first. 



3 Adiabatic approximation 

In this section we summarise the main points of the adiabatic approach and establish their 
connection to some facts about the parametrically driven, damped NLS equation, available 
in literature. Consider equation ([!]) with TZ as in (Q), and assume that h and 7 are small. In 



the adiabatic approximation the soliton solution is assumed to be given simply by the unper- 
turbed NLS soliton fl2lD with the parameters £, r], a, and 6 being slowly changing functions 
of time. We now derive and discuss the adiabatic equations for the soliton parameters. 

The matrix T(£) = T(— oo, oo; (), the key element of the perturbation theory defined in 
(|23|) , can be easily computed for the NLS soliton with the t-dependent parameters. In this 
case = I and *+(C) = T(C). Using ([18]) we get 



Y 22(C) = 7 / d^sech z 



e e 
+ 



c - Ci c - cr 



(30) 



and 



T«(C) 



dz exp(— 2iCa; + i<y?)sech 2 2; 



(C-CDe 22 (C-Ci)e- 22 1 
2i V (( - Ci) 2i V (( - CD i77 



2ir)Kl(z) 



+ (c-Ci)(c-cr)J- (3 j 

In equations fl30|) and (|31|) we have omitted, for notational convenience, the explicit t- 
dependence and defined 

n = e- ilp TZ. (32) 

The variables z and <p are defined as in (|T9| ) where this time, we need to take into account 
the t-dependence of the zero Ci = £ + ^7- I* 1 this case the integration of (|T6|) gives equation 
dTTj) with 

h = K x x-2i [ dtd (33) 
whence the core phase and position of the soliton are found to be 



ifo = e + 2£x + 4 J ' dt(rf - e 2 ) , x = ^ a + 8 J dt^rj 

V 

Inserting (|30| ) and fl3~T|) into (|26| ) and ( p7|) (with j = 1), gives 

00 

+ = _i / dze z sech 2 z[ft (z)+7^(-z) 



(34) 



(35) 



and 

A 

dt 



[a+it 



a+8 / dt^rj 



dt 4?7 



dz sech 2 z ^cosh z— ze* 



, (36) 



where, as in the previous section, (i = £ + it], and we have used the identity 



-2idx + icp = i6 + a-z-Ai / dt 



(This follows from (0) and @.) 

Let TZq s stand for the value of the quantity TZo calculated on the soliton, i.e., with the 
phase and position parameters given by (Rl). Equations fl3"o] ) and with 7Z = lZ 0s 
constitute the set of the adiabatic equations for the soliton parameters. These equations are 
equivalent to those derived by Karpman and Maslov EH . 

For the parametrically driven, damped NLS equation (|])-(@) the quantity 7Z 0s equals 



e- lv {he 2inot q* s 



vyq t 



2?7 1 7 — ih exp 2i(Q t — f) | sech^;. 



Substituting this into ( |35l) and (|36|) we obtain, after some algebra: 

d.77 2nh£ 



dt smh(7i rj) 
d£ 2nhe/v 



d / 

— la + it 
dt 



dt sinh(7r£/?7) 



sin(2fi t — 2(^ ) — 2777, 
sin(2fio^ — 2y?o); 



(37) 

(38) 
(39) 



sinh(7r£/?]) 



cos(2fi /; — 2(y9 ) 



— n cothf— 



z 

V 



+ ^x — (ri-i^). 



(40) 



Equations (|38l ) and ([39]) imply d(?7£)/df = —2^1]^, i.e the quantity 7?£, proportional to the 
momentum of the soliton pi| , has to decay to zero as t — > 00. Taking this into account we 
will restrict ourselves to the nonpropagating soliton: £ = 0. We can also choose a = so 
that the soliton is placed at the origin: x = 0, see equation (]33J). Thus, sending in ( |3lf ) and 
(|40D £ — ► and making use of the identity 



lim ^ N 



111 

— 77 COth 

S V 77 



(f)0- 



we arrive at a closed system of equations for the soliton's amplitude and phase. This can 
be presented in the following convenient form: 



i) = — 277(7 + h sin $), 
$ = 2fi - 8?7 2 - 2/icos$, 



(41) 
(42) 



where 



$/2 = ft i - (43) 
is the difference between the phase of the driver and the core phase of the soliton. 



These equations were first derived in |36j and [J57I] (within a different approach, though). 
The two first-order equations ( |41]) -( [42f) can be rewritten as a single second-order equation 
for $: 

$ + 2^27 + /isin$^6-8^O -/icos$^7 + /isin^ =0, (44) 

which in some situations is more amenable for analysis. Equations (f4~T]) -(42) have two fixed 
points (j]±, $±) which correspond to stationary soliton solutions ( pT| ) with i] = r]± and 
$ = $ ± . Here 

2?7 ± = \/fl ±H, $ + = -7r + arcsin , $_ = - arcsin (45) 

and 

H = v//i 2 - 7 2 . 

Although obtained in the adiabatic approximation, these two solitons turn out to be exact 
solutions of the parametrically driven damped NLS |T3|] . 

Linearising equation (|44|) in the small perturbation <p = $ — $±, and making use of 
relations 

/isin($ ± + 0) = - 7T #0 + £>(0 2 ), /icos($± + 0) = T# + 70 + C(0 2 ), (46) 
yields the equation of damped linear oscillator: 

4> + 2 7 ± 8H(n + H^(p = 0. 

From here one readily concludes [37]] that the g + -soliton (i.e. the soliton (pl|) with r] = r] + 



and $ = $+) is adiabatically stable and, when excited, exhibits decaying oscillations at the 
bare frequency fi s = a/8//(^o + H). The g_-soliton (for which r\ = r\- and $ = $_) is 
adiabatically unstable, and this of course implies that it is unstable within the full partial 
differential equation. (This is corroborated by the stability analysis of the full equation, see 
||13|| .) For this reason we disregard the g_ soliton in what follows and focus entirely on the 
q + ; from now on the "parametrically driven damped soliton" will always mean the soliton 

q+- 



It is interesting to note that the adiabatic equations (^)-(^) have another solution that 
admits a simple interpretation. It is given by rj = and $ = 2(Qq — h cos $) and corresponds 
to the flat solution q = of the full damped-driven NLS. It will reappear in our analysis of 
the soliton-radiation interaction below. 

The main shortcoming of the adiabatic approximation is in that it ignores the soliton- 
radiation interaction. As a result, it is unable to reproduce many features of soliton's 
dynamics even for fairly small perturbations. In particular, the adiabatic approach does 
not capture the oscillatory instability of the q + soliton which sets in as the driving strength 
exceeds a certain — rather low — threshold (For example, for 7 = this threshold 



is at h = h c (0) = 0.064.) Neither is it capable of reproducing secondary bifurcations and 
chaotic dynamics of the soliton. In what follows we go beyond the adiabatic approximation 
and take the soliton-radiation coupling into account. 



4 Evolution of the spectral data 



4.1 The essence of our approach 

As we mentioned in section 2, the soliton of the unperturbed, "pure", NLS corresponds to 
a single zero of the Riemann-Hilbert problem associated with this equation. Our approach 
to the parametricaly driven, damped NLS is based on the fact that it can be cast in 

the form 

iu t + u xx + 2\u\ 2 u = he 2int u* - (H + ij)u, (47) 

where 

u(x,t) = e iHt q(x,t), (48) 

and Q = Qq + H. The key property of the new formulation is that both the left and the 
right-hand side of equation (|47|) are equal to zero for u equal to 



u+(x,t) = 2irj + != — r (49) 

cosh {27] + x) 

where $ + and r/ + are as in (|45|) . (In particular, 2rj + = y/Q.) Hence u + coincides — exactly 
- with the soliton of the unperturbed NLS equation with a particular amplitude and phase 
(selected by the parameters of the perturbation). Therefore, this solution of the perturbed 
equation is also associated with a single zero of the Riemann-Hilbert problem underlying 
the unperturbed, integrable, NLS. It is important to emphasise that this correspondence is 
valid for arbitrarily large values of h and 7 (where h only has to be greater than 7 so that 
the soliton flU), (H) exists). 

Rescaling t,x,q,h and 7 we can always arrange that Q = 1 in equation (||). Next, we 
have already seen that the motionless soliton is a solution of the adiabatic equations. It is 
not difficult to realise that, in our case, even when the radiations are taken into account, 
an initially quiescent soliton will remain nonpropagating at all times. This follows from 
the evolution equations for the spectral data {r/, £, xo, <po] g(C)} for the damped-driven NLS 
([[])- (|2|). These evolution equations, linearised in (p and g((), can be easily shown to be 
compatible with the constraint 

x = 0, £ = 0, g(-Q=g(Q. (50) 

Consequently, in this paper we confine ourselves to the internal dynamics of the nonprop- 
agating soliton (placed at the origin for convenience) and its radiation. The corresponding 
solution of equation ( P7| ) is given by an even function of x: 

u(x,t) = u s (x,t) + u r (x,t), (51) 

where u s (x, t) has the form of the motionless soliton of the "pure" NLS, located at the origin, 
with the time-dependent amplitude and phase: 

u a (x,t) = 2w?e l ^°sech2. (52) 



Here 77 = rj(t), z = —2r)(t)x, and ipo is given by equation ( |34] ) where we only need to set 
^ = 0: 



Equivalently, the soliton can be written in the form 

. . exp {ifit - i$/2} 

u s {x,t) = 2l7] , 

cosh ,2 

where the variable $ is defined by 

$/2 = ttt - <p . (53) 

The definition (|53|) is equivalent to (0); in both cases $/2 is the difference between the 
phase of the driver and the core phase of the soliton. (Here we should alert the reader to 
the fact that, since q and u are different by the factor (f48|), the core phases <fo(t) in (|53|) 
and in (^) are different by Ht.) 

The second term in (|5T| ) accounts for the radiation waves. As in Ref. P3| , in our deriva- 
tion of the evolution equations for the spectral data we will retain only terms linear in 
radiation. Hence it is sufficient to solve the linearised version of the regular Riemann- 



Hilbert problem (|15|) to obtain the radiation part of the solution (0). The linearisation of 
the Riemann- HUbert problem produces the Plemelj jump problem: 

_ = 1V'-'"'Y,V '-'"T ; . ImC = 0. (54) 

Taking into account the normalisation condition ty ~~ > I as C ~~ > 00 ) we obtain from 



= 1 + h I rh m ( c2 "o 9(,) ) r " (,) - (55) 

—00 

(Here the sign + respectively — indicates that ( lies in the upper respectively lower half- 
plane.) The radiation part of the solution (pT|) is now given by equation (El): 



00 

u r = -2 lim C(^o + )i2 = - / dC [ruCr- 1 )^ e 2l ^g + T 12 (r _1 ) 12 e" 2i V] . (56) 

—00 

Finally, substituting the one-soliton dressing matrix (|I8|) into equation (|56"D, introducing the 
notation 

fc = - (57) 
77 

and defining the radiation amplitude 6(/c) by 

b(k) = exp(-i<p )g(r}k), (58) 

we arrive at the formula for the linear radiation: 

00 

uJx, t) = 1— \ -f — Uk - i tanhz) 2 e~ ikz b(k) + sech 2 ze ikz b*(k)] (59) 

TTl J fc + 1 L J 



in exact agreement with the corresponding result in Ref. [25 



4.2 A closed system for the evolution of the spectral data 



Since the r.h.s. of the perturbed NLS (|4"T| ) is linear in u, the decomposition of the solution 
into the soliton and radiation parts, equation (|51~D , induces the corresponding decomposition 
of the perturbation: 

n = tz(u s ) + n{ur) = n s + n r . 

The perturbation matrix fl22| ) splits accordingly: R = R s + R r . Substituting for u s and u r 
from fl52"|) and (|59|), we get 



TZ S = (R s ) l2 = he 2int u* s -(H + ij)u s = (-iH + 7 - ihe i<s> ) 2r]e itpo sechz 



and 



TZ r = (R r ) 12 = he 2im u* - (H + ij)u r , 

respectively. Discarding terms higher than linear in b in the expansion of the functional (^) 
yields 

00 

T(*) = / dxexp(-*^ 3 )r-(^){i i , + « r + [fl.,*„ + (^)]}r^)ex P (*^3), (60) 



where T(() is the one-soliton dressing matrix (|TJ) and ^+(0 is given by equation (|55|). 

In order to obtain the evolution equations for the spectral data one has to evaluate the 
integrals in the expression (|60|); substitute the resulting matrix T in equations d26|)-(p8|) 
and use the simplifying conditions (|50D for the nonpropagating soliton. (As we already 
mentioned, these conditions are compatible with the evolution equations for the spectral 
data.) It is convenient to present the final result in terms of $; w = 8rj 2 ; and the real and 
imaginary parts of the radiation amplitude b(k) = 5r + ibj. After tedious but otherwise 
straightforward calculations one arrives at the following system: 



w = — < (7 + h sin $) 



4 — / dfcsech 



$ = 2(1 — /icos$) - w, 



(61) 



b R (k) 



— 2h cos $ / d/csech 



71k 



h{k) \w 



db 



R 



ot 



w 



-(l + A; 2 ) + 2/icos$(/ + /C) 
2 



bi - (7 + h sin $) 



Ob 



(62) 



b R + 2k^+[M+- M\b 



dk 



— (7 + h sin $)7rsech 



^ = -|(1 + k 2 )b R - 276/ + (7 + h sin $) 



bi - 2k 



dbj 
~dk 



M-)bj 



(63) 
(64) 



Here the operators /C and Xi± are defined on even functions: 

(l-k) (I -Ik) 



£/(*) 



d/ 



sinh[vr(/ - fc)/2] (1 + F)(1 + / 2 



«m 1 ! a, ( l±k ) (l 2 + k 2 + 2) 
• M±/(fc) = 2 J <U 8 tohKl-fc)/2](l + f)(l + P) / " ) - 

— OO 

where the singular integral in the expression for .M+ should be understood in the sense of 
the Cauchy principal value. 

In equations ( JjHD and (0), the notation db/dt is meant to indicate that the derivatives 
are taken for fixed k. (On the contrary, writing b would mean db/dt + kdb/dk with k(t) as 
in (|57p, A; = £/rj(t).) We are using partials here for later computational convenience. 

It is worth re-emphasising that equations (|6TD-(|64|) are valid for arbitrarily large h and 
7. The only approximation we made in their derivation, was to drop terms of order higher 
than linear in b. 

Letting b R = bj = reduces equations (|6lD- (|62D to the adiabatic equations ([iT|)- (|42]) . 
Like the adiabatic equations, the system (]6"T|)-(|6"4|) has a fixed point w = 2Q = 2(1 + H), 
$ = $ + = —7i + arcsin(7//i), bn = b I = 0, which corresponds to the soliton ([49]) of equation 
(f4"Tp. Another meaningful solution arises by setting w = and solving floTp for $(t); the b^ 



and 6/ are then recovered from the nonhomogeneous linear system fl63|)-064|). We will return 
to (a descendant of) this solution in section 5. 

4.3 Linearised evolution of the spectral data 

Linearising the system (|61~D-(|64D about the above fixed point, one obtains four first-order 
linear equations which are equivalent to a pair of second-order equations for <fi = $ — $ + 
and P(k), where 

m = ^^ttL (65) 
The second-order system has the form: 

00 

<j> + 2j<j> + SHQcf) = 2 V^HQ J dkVl + k 2 sech (3(k), 

—00 

(66) 

p + 27/3 + Q/3 = -2^HQVl + k 2 sech J 0. 

Here we have made use of relations fl46|) and introduced a real symmetric operator Q: 

Q = (1 + fiA; 2 ) 2 - # 2 + 2iy^B, (67) 

where 

00 

J sinh[jr(fc - /)/2] V(l + * 2 )(l + ( 2 ) 

— OO 

It is worth noting here that since the equation for is not uncoupled from the equation for 
(3, the 4> and /3(k) do not represent the normal modes of the soliton-radiation system. This 



is of course a consequence of the nonintegrability of the perturbed nonlinear Schrodinger 
equation. 

The system (|fi?S|) is exactly equivalent to the NLS Q37D linearised about the soliton P5|). 
In particular, one can readily check that for H —* 0, the oscillation frequencies of (|66|) 
reproduce the asymptotic expansions obtained in |I3||. Indeed, letting <p{t) = e^ -7 ^ and 
P(k; t) = ^ iu) -^ l a(k) transforms fl6|) to an eigenvalue problem 



oo 

iHQx - ZV^Htt J dk Vl + k 2 sech ( y J 



(69) 



2v / vr#Wl + /c 2 sech 
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where [l = oo 2 + r ) 2 . (Below, we normalise the eigenvector {x, a(k)} by setting \ - 
A simple perturbation analysis shows that for small H there is an eigenvalue [i 



l.: 



8#[1 + 



(7r /6 — 1)//] + 0(H ). Another discrete eigenvalue detaches from the boundary of the 
continuous spectrum (which extends from fx = 1 — H 2 to infinity): fi = 1 — 9 if 2 + 0(H 3 ). 
Both eigenvalues coincide, to within terms 0(H 3 ), with the corresponding expressions for 
the eigenfrequencies in the space of fields u(x,t) |13| . 

In fact the eigenvalue problem (|6~9D is amenable to analytical study not only in the H — > 
limit. (This is the principal advantage of the analysis in the space of scattering data over 
the linearisation in the w(:r, t) -space where one has to resort to the help of computer.) The 
radiation can be diagonalised in the basis of eigenfunctions of the operator Q: 



Qip{k) = Ei){k). 

This is equivalent to the following differential eigenvalue problem: 



-H 2 -E 



y = AHVt 



( sech 2 £ — ^ ) + sech 2 £ y 



(70) 



(71) 



Here y = ?/(£) is the Fourier cosine transform of the function (1 + k 2 ) 1 ^ 2 ip(k): 

oo 

y(0= [ 4= «(*e)-^- P2) 



'2ti 



VTTk 2 ' 



The continuous spectrum of the operator Q occupies the semiaxis E > 1 — H 2 while discrete 
eigenvalues satisfy E n < 1—H 2 . A more accurate bound on discrete eigenvalues is established 
in the Appendix: 

^5 HVL 



'\-H l 



21- # 2 



<E n <l-H z . 



(73) 



Equation ( |73|) implies that the operator Q cannot have any discrete eigenvalues for H = 0. 
However, there is at least one discrete eigenvalue for any positive H (see the Appendix.) 

Expanding the radiation amplitudes a(k) over the orthonormalised basis of (even) eigen- 
functions of Q: 

oo 

a(k) = J2cn*p En (k)+ J dEc(E)ip E (k), 
the eigenvalue problem (^9|) is cast in the form 

oo 

(8#0-/i) = J2< E n)c n + J <1Ek{E)c{E), (74) 

(fi — E n )c n = K,(E n ), (75) 
(ji-Bf)c{E) = K(E), (76) 
where k measures the coupling of the soliton and radiation modes: 

oo 

k(E) = 2y^HQ j dk VTT¥ sech H^-j ip E (k) 

— oo 

oo 

= 2V2Hil f deseche^l-^ y E (0- 



(77) 



Here Qip E {k) = E%fj E {k)] ip(k) and ?/(£) are related by the Fourier transform (|72|). 

Solving (|75|) and (|7"6[) for c n and c(E), respectively, and substituting these into ([74]) we 
arrive at the characteristic equation of the form 

fi-BHQ = s{p), (78) 

where the function 



oo 

M=9^,H) = J2P^+ [ dE 

E n - fi J 



k 2 (E) , . 



1-H 2 

is defined for fi < 1 — H 2 . This is a monotonically growing function for all \i except points 
fj, = E n where it drops from plus to minus infinity. For \i < E% the function is strictly 
positive and decays to zero as \i — ► — oo. (See Fig.p. Roots of equation ( [78] ) give discrete 
frequencies of the system (|66|), while the spectrum of continuous frequencies extends from 
jj, = 1 — H 2 = 1 + 7 2 — h 2 , to infinity. 



o 



1-H 2 



Figure 1: The graphical solution of the characteristic equation ([781). The solid straight line 
depicts the function y(/j) = fi — 8HQ and the two solid curves intersecting it are branches of 

y = sO)- 



Before discussing the roots, we would like to note that the characteristic equation (ff8|) 
does not contain 7 explicitly. Having solved it for fi we can recover uj for any 7: uj = 
\ffa- 7 2 - This invariance of the eigenvalue problem in the space of scattering data is 
an exact equivalent of the invariance of the linearised eigenvalue problem in the space of 



solutions to the NLS equation [|13[ . 

The graph Fig.|TJ shows that the function y = g(/i) and the straight line y = fa-8HQ have 
at least three intersections, fa, fa and fa such that fa < fa < E± < fa. As H grows, the 
straight line is shifted down and the function g(/i, H) also changes. In particular, for small 
H the coefficients k 2 increase in proportion to H 2 while the eigenvalue E\ decreases and the 
edge of the continuous spectrum, 1 — H 2 , also moves to the left. Hence for any fixed fi < E\ 
the corresponding grows. As a result, the intersection points \i\ and approach each 
other, then merge for some critical value H = H CT and emerge into the complex plane. As 
H is increased past H CT , the imaginary parts of fa = /4 grow and eventually the positive 
imaginary part becomes equal to 7. This is the point of the Hopf bifurcation; for H above 
this point the soliton (^) is unstable. 

Conclusions of the above graphical analysis are in exact agreement with the behaviour of 
eigenvalues of the linearisation of equation (f|7D observed numerically |13| . A new feature is 



the appearance of another discrete eigenvalue, fa, which was not detected in the numerical 



computations of [13|, |18| . One reason for this omission could be that the separation of fa 
from the continuum remains exponentially small for small H. Indeed, for H — > equation 
(|77|) infers that k^) = AirH 3 / 2 + 0(H 5 / 2 ) and k(1 - H 2 ) = V^H + 0(H 2 ). Assuming 
fa = l — H 2 — 5fa where dfi/H 2 — > as H — >• 0, equation (f79|) gives g(/i, H) = —8nH 2 \n5fa- 



2n 2 H + 0(H 2 ). Substituting into flT8|), we obtain for the separation 5p: 

^ = exp {-^-^i + c } as/ ^°- (80) 

where C = C(H) remains bounded as H — > 0. As H grows, the frequency /J3 remains real 
and therefore does not give rise to any instabilities of the stationary solitons. However it 
may play a role in the resonance phenomena involving oscillating solitons. 

We conclude this section with two remarks. First, our linearisation in the space of 
scattering data allows to explain the origin of the oscillatory instability of the soliton. When 
H is very small, the soliton oscillations are virtually uncoupled from the radiation waves. 
The former have their frequency close to \i — while the radiations have continuum of 
frequencies occupying the semiaxis /x > 1. As if is increased, the coupling grows, and as 
a result of that the soliton's frequency is dragged closer to the continuum while another 
local mode is pulled out of the radiation spectrum. Finally, the two modes merge and the 
resonance occurs. 

Second, we observe that as k grows, the right-hand side of the bottom equation in (|66|) 
decreases rapidly. We also notice that only /3(k) with small k contribute significantly to 
the right-hand side of the top equation. Consequently, only long radiation waves couple to 
the soliton. In the next section we study the effect of the k = radiation on the soliton's 
nonlinear dynamics. 



5 The long-wavelength limit 

In order to study the soliton-radiation interaction at the long- wavelength limit, we write 
the spectral density b(k) in the form b(k) = b e (k) = e~ 1 b{k/e), where b is an even function 
of its argument; integrate equations (|63|) and (0) over a small range of k G [—5k, 5k] and 
then send e — > in (pT|)-(p^). This allows to derive a finite-dimensional system without 
any ad-hoc cut-off parameters (cf. [^, [34]]). The resulting four- dimensional system 
comprises equations for the amplitude and phase of the soliton, and the complex amplitude 
of the k = radiation: 

8k 



p + ia = lim J dkb e (k). 



-6k 



Skipping some lengthy but elementary calculations, we produce only the final result: 

$ = 2(1 - /icos$) - w, (81) 

w = -{(7 + /isin$)(4 - p) - 2a/icos$|w, (82) 

p= |^ - cH+ (c + 2)/icos$ja- (7 + /i sin $)(4 - p), (83) 

o = (~ + ^ + i/icos$ ) p+ (7 + 3h sin $)cr. (84) 



Here w = 8r) 2 and c = 41n2-2« 0.77. 

Below we will compare conclusions based on the analysis of the four-dimensional system 
flgTD-flSip with solutions of the full nonreduced NLS equation. To facilitate the comparison 
with results available in literature, we take the damped-driven NLS in autonomous form: 

#t + ^ + 2|#ty = + (85) 

Here if) is related to a solution of equation ( f4T|) by a simple phase transformation: ip(x, t) = 
— exp(— iQt)u(x, t). The value of the field (which is a sum of the soliton and radiation) at 
the centre of the soliton, i.e., at x = 0, can be expressed through the variables of the reduced 
system. Using (|5T|) and (|59"D we obtain in the limit e — ■> 0: 



u(x, t) 



x=0 

and hence the value of the solution of the NLS (EB) at the point x = is given by 



ip{x,t) 



x=0 

The finite-dimensional system has two exact solutions. One of these, 

$ = -7r + arcsin(7// i ), w = 2(l + H), p = a = 0, (87) 

corresponds to the pure-soliton solution of equation (^). The second solution arises by 
letting w = 0, with $ determined from the equation (0), and p and a defined by the linear 
system fl83|) - (0). Choosing the constant of integration so that $(0) = 0, we obtain 



$(t) = 2arctan <^ tan(Vl - h 2 t) } mod (2tt) 




Substituting this into equations fl8"3|) and (|S4"D, the system for p and cr becomes 

p = ( 7 + /id 1 )( / o-4) + (41n2/id 2 -c)(T, 



59) 



cr = (7 + ?>hdi)a + ~(H - hd 2 )p, 



where 



VI - h 2 sin(2Vl - fc 2 *) , ^/x /» + cos(2 v / T 3 7?t) 

tZi = sin $(t) = v ; , d 2 = cos $(t) 



hcos ' 1 l + / i cos(2 v / r^7?t) 



In view of equations (|52| ) and (|59"D, u> = amounts to ijj(x,t) = u(x,t) = and hence 
the above w = solution corresponds to the flat zero solution of the full nonreduced NLS 
equation (p5|). This fact is quite remarkable. Indeed, although the system (|8~1"|)-(|84|) was 



derived under the assumption of the proximity to the (damped driven) soliton, it does 




Figure 2: A numerical evolution of w, p, and a: an attraction to the w = solution and an 
explosive instability of the latter due to the exponential growth of p and a. 



possess the flat solution which cannot be regarded as the soliton's small perturbation. As 
in the full partial differential equation (j85|), where the unstable soliton can decay to the flat 
attractor, finite-dimensional trajectories starting near the unstable fixed point (^) can be 
attracted to the w = solution. An example of such evolution, obtained numerically, is 
presented in Fig. |, for t G [0, 250]. 

It is important to note here that for 7 > 0, solutions to the linear system (^) are 
exponentially growing functions (see Fig. @(b)). This is not in disagreement with the fact 
that the corresponding ip(x,t) is zero for all x and t. Indeed, the total amplitude of the 
long-wavelength radiation includes a factor 77 (see e.g. ([5£]) or ([36l).) Therefore, if w = r\ = 0, 
the total amplitude is zero no matter what p and a are equal to. 

The fact that p and a are exponentially growing functions, gives rise to an instability of 
the w = solution of the system fl%l"|)-(pi|). Indeed, the multiplier 

x(t) = - {(7 + /isin$)(4 - p) - 20-/icos$} 

in the right-hand side of (|82"D , may assume large positive values. (For example, each time 
sin$ goes through 1, x{t) becomes equal to approximately (7 + h)p.) This means that for 
w^l we have two different time scales in the system: the variables $, p and a do not 
change appreciably on intervals At < 1 and can be considered constant, while w grows with 
the exponential growth rate t ^> 1. Consequently, if w is assigned a small but nonzero 
value at the moment of time when x(t) is large, it will quickly (within At ~ t _1 ) grow to 
values of order x. This is indeed observed in simulations, see Fig.^| for t > 250. However, 
the instability of the w = solution is spurious — in the sense that it does not mirror any 
genuine instabilities of the zero solution in the full NLS equation. 

6 Reduced finite-dimensional dynamics 
6.1 Onset of instability 

Linearising the system (|81|)-(|84l) about the fixed point (|87D and using (|46|), we obtain 

<j> + 2 7 + 8H(1 + H)(f) = AH(1 + H)a, (90) 

a + 27a + (1 + H)\l - (2c + l)H]a = -4#(1 + H)<p. (91) 

Here (f) is the perturbation of the soliton's phase, defined as </> = $ — $ + . Letting <p(t) = xe xt 
and a(t) = o"oe At , where A = —7 + iui, we obtain the characteristic equation for complex uj: 

u) 4 - (tt 2 s + Q 2 r - 2 7 V 2 + {Q 2 S - 7 2 )(fi 2 - 7 2 ) + 16# 2 (1 + Hf = 0. (92) 

Here 

ft 2 = 8H(1 + H), n 2 r = (l + H)[l- (2c + l)H\. 

Like the eigenvalue problem (|7^ ) - ([76|) for the full set of spectral data (and like the 
eigenvalue problem for the underlying partial differential equation [0]), the characteristic 



7 2 + uj 2 : 



equation (|92| ) can be conveniently reformulated in terms of the self-similar variable /i = 

fi 2 - P fi + q = 0, (93) 



p = p(H) =Vl 2 s + ft 2 , q = q(H) = ft 2 ft 2 + 16H 2^ + R y 



where 

The fact that the reduced system (|81~|)-(|34"|) inherits the self-similarity of the parent PDE, 
deserves to be specially emphasised. It implies that the reduction procedure based on 
the Riemann-Hilbert problem — unlike the variational reductions |42], |43| - preserves the 
structure of the infinite-dimensional phase space. 

The fixed point is unstable when ReA > 0. In terms of /i, this condition translates into 

[lm^(H)} 2 ^ , 



4Re//(#) 



> 7 • (94) 



The discriminant of (OBf) is 



D(H) = (HiH 2 ) (1 + H) 2 (H - Hi)(H - H 



2)- 



where E x = (8 In 2 + 13)" 1 w 0.054 and H 2 = (8 In 2 - 3)" 1 w 0.39. Between H t and H 2 
the quadratic (p^) has a pair of complex-conjugate roots while outside this interval both 
roots are real and nonnegative (Fig. 3). Consequently, the inequality ( j94l) can only hold for 
Hi < H < H 2 . For H within this interval the inequality (0) amounts to 

8 7 2 < -D(H)/p(H). (95) 

For small 7 the cubic equation D(H)/p(H) = — 87 2 has one negative and two positive roots 
which we denote Hi(y) and H 2 (y), < Hi (7) < H 2 ( y). (No te that 7^(0) = H x and 
7i 2 (0) = H 2 ; hence the notation.) Recalling that H = yjh 2 — 7 2 , the instability inequality 
(|95D can be rewritten as 

^?(7)+7 2 <^< V / ^i(7)+7 2 - (96) 

As 7 is increased, the positive roots ?ii(y) and Ti. 2 {y) merge and become complex. For 
greater 7 the inequality (|9~4D-(|9~5D cannot be satisfied by any h. 

Equation (|96| ) gives an explicit form of the instability region on the (h, 7)-plane (Fig.4). 
As h is increased past the value h = \/Hl + 7 2 (or decreased below h = ^H 2 + 7 2 ), the 
fixed point becomes unstable via a Hopf bifurcation. 

The instability setting in for H > ^1(7) corresponds to the Hopf instability of the 
soliton in the full damped-driven NLS equation. Note however that the finite-dimensional 
instability threshold h = \fH\(y) + 7 2 is somewhat lower than the instability threshold in 
the partial differential equation (which is also shown in Fig. 4). On the other hand, the 
upper boundary of the instability region, h = \jH 2 {y) + 7 2 , does not have a counterpart in 
the full NLS equation. (In the full PDE, the soliton does not restabilise as h is increased 
This fact alone is sufficient to conclude that the finite-dimensional system cannot be 




Figure 3: The real (a) and imaginary (b) parts of the roots to the characteristic equation ( |93|) . 




Figure 4: Solid curve: the boundary of the instability region of the fixed point. (The fixed 
point is unstable inside the "parabola".) Dashed curve: the line of the Hopf bifurcation of 
the stationary soliton of the full NLS equation (|85|). (The soliton is unstable above the dashed 
curve.) Also shown is the straight line h = 7, the boundary of the existence domain of both the 
soliton and the fixed point ( |87|) . 



expected to provide a good approximation to the infinite-dimensional dynamics for h greater 
than approximately 0.4. 

In order to find attractors in the region inside the parabola (|96|), where the fixed point 
is unstable, we performed a series of numerical simulations of equations fl81"D-(84]). Here our 



strategy was similar to that used in the simulations of the full nonreduced NLS equation [|17 
that is, we varied h for a fixed value of 7. We also adopted the same strategy with regard to 
the choice of the initial conditions. Our simulations always started with the (unstable) fixed 
point perturbed by values of order 10~ 14 which is several orders of magnitude smaller than 
the local discretisation error of the Runge-Kutta approximation. (For some values of h and 
7 the final state of the system was extremely sensitive to tiny changes of this perturbation.) 



6.2 Finite-dimensional attractors; 7 7^ 

As h is increased for a fixed nonzero 7 (with 7 < 0.31), a limit cycle is born supercritically 
at the point of the Hopf bifurcation, h = yR^y) + 7 2 - This is in exact agreement with the 
full partial differential equation. The subsequent bifurcation diagram depends on the value 
of 7. 

For 7 greater than approximately 0.2, the stable limit cycle persists as h is increased all 
the way up to h = a/^Kt) + 7 2 where it shrinks back to the fixed point. On the contrary, if 
we increase h for a smaller 7, 0.12 < 7 < 0.20, the limit cycle looses its stability at a certain 
value h = h z (where h z < \J1~L\ + 7 2 .) For h above h ZJ the finite-dimensional trajectory 



emanating from any initial condition, quickly settles to the w = solution corresponding 
to the zero attractor of the full nonreduced NLS. (This is illustrated by FigfJ and, for a 
different 7, Fig.^(b).) Increasing h still further, the w = solution persists as the only 
attractor over a sizeable range of h values — until a stable limit cycle reappears. This range 
becomes wider for smaller values of 7. 

The range of driving strengths where the only attractor is w = 0, exists for all 7 G 
(0, 0.20] (although for some 7 the limit cycle may undergo a number of intermediate bifur- 
cations before disappearing from the attractor chart; see below.) This is in exact agreement 
with the behaviour observed in the full partial differential equation JET] where the solitonic 
attractor undergoes a crisis and the flat zero solution remains the only attractor in the sys- 
tem. Thus we may conclude that taking into account the coupling of the soliton to the k = 
radiation, is sufficient to explain the occurrence of the "desert region" on the (h, 7)-plane. 
It is appropriate to mention here that a similar "desert" spanned just by a flat attractor, 
arises in the externally driven NLS El 133. It is natural to assume that the appearance of 



the "flat desert" in the latter system can be still explained by the soliton-longwave radiation 
coupling. (The failure of the four- dimensional reduced systems proposed in [^, [38], |34| to 
capture the crisis of the localised attractor and reproduce the "desert" , should be probably 
attributed to nonoptimal variational Ansatze.) 

Although our finite-dimensional system with 7 < 0.2 correctly reproduces the sequence of 
attractors arising in the full PDE (stationary soliton — > oscillating soliton — > zero solution), 
it does not necessarily capture fine details of this sequence. Unlike the oscillating soliton in 
the full NLS equation, for 7 between 0.12 and 0.20 the finite-dimensional limit cycle does 
not undergo any period-doubling bifurcations. The largest 7 for which the period-doubling 
occurs in the reduced system (^1|)-(|^), is 7 = 0.12. In this case the period-2 cycle arises 
at h = 0.16 and then degenerates back to the period-1 as h is increased past h = 0.17. 
Increasing h still further, the period-1 yields to the w = solution for h > 0.195, without 
any intermediate period-doublings. A similar pattern arises for 7 = 0.11. 

In the case 7 = 0.10 the sequence of finite-dimensional attractors is richer. In this case 
we observed the whole cascade of period-doubling bifurcations, culminating, for h = 0.146, 
in a chaotic attractor (centred on the fixed point). For h between 0.146 and 0.150 the only 
attractor was found to be w = 0; for h = 0.150 the chaotic attractor reappears, and for 
even greater h it degenerates to the period-2 (for h = 0.165) and then period-1 cycle (for 
h = 0.170). At h = 0.18 the cycle disappears and w = was the only attractor we detected 
in a wide band of h values above h = 0.18. Next, increasing h for the fixed 7 = 0.05 and 
0.07, the period-1 (PI) limit cycle yields to the period-2 (P2) and then to w = 0. For 
7 = 0.06 the sequence of attractors was PI — > P2 — > (w = 0) — > P3 — > P4 — > (a 4-band 
chaotic attractor) — > {w = 0) — » (a 2-band chaos coexisting with P3) — > (a 1-band chaos) — > 
(w = 0). (The last two regimes are illustrated in Fig.^.) The full NLS equation also exhibits 
a rich host of attractors for 7 > 0.05, including the chaotic soliton 1 17] , [3D]; however, details 
of the two bifurcation diagrams do not necessarily coincide. 

It is only for very small 7 that the finite- and infinite-dimensional dynamics are in exact 
agreement. Increasing h for the fixed 7 = 0.01, 0.02 and 0.04, the period-1 cycle yields 
directly to the w = attractor. This is consistent with the absence of the period-doubling 
in the full PDE. The smallest 7 for which the period-2 soliton was observed there, was 
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Figure 5: Trajectories in the four-dimensional system (pl|)-(p4D with 7 = 0.06. (a) A strange 
attractor corresponding to the chaotic soliton. (b) A crisis of the chaotic attractor: the trajectory 
escapes from the neighbourhood of the fixed point and flows into the w = solution. To facilitate 
the comparison with the corresponding regimes in the PDE, the trajectories are presented in terms 
of the real and imaginary parts of the quantity 



7 



0.055 [01. 



6.3 The undamped case, 7 = 

This case deserves a separate consideration as the arising attractors are very different from 
those occurring for nonzero damping. Consistently with the PDE [JTjJ, no limit cycles were 
detected for 7 = 0. Instead, we observed two types of chaotic regimes. For h in a narrow 
interval 0.0539 < h < 0.06 just above the Hopf bifurcation value h = Hi = 0.0539, the 
trajectory was seen to wind chaotically in an annulus centred on the unstable fixed point 
(Fig.|6|a). The outer radius of the annulus is approximately one order of magnitude smaller 
than w = 2(l+h) , the w-coordinate of the fixed point. This proximity of the orbit to the fixed 
point justifies one of the assumptions made in the derivation of the finite-dimensional system 
(|8T]) - (p4|) . The strange attractor shown in Fig.|6|a represents small chaotic oscillations of the 



amplitude of the soliton about its stationary value l^l^o = \A + h. The chaotic oscillations 
of the undamped soliton were indeed observed in numerical simulations of the full nonreduced 
NLS equation |40| ; however these would die out as transients and the evolution settle to 



either a slowly growing or decaying breather [ITS). The fact that the chaotic attractor does 



not arise in the full partial differential equation with 7 = and only persists as a transient 
chaos, can be explained by the interaction of the soliton with short and medium wavelength 
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Figure 6: Trajectories in the finite-dimensional system (|8lD-(|84D with 7 = 0. The strange 
attractor (a) and its crisis (b). 

radiations — which we have neglected in our present derivation of the finite-dimensional 
system. For h greater than 0.06 (but smaller than H 2 = 0.39, the upper boundary of the 
fixed point's instability domain), the chaotic solution is no longer centred on the fixed point. 
(See Fig.^D). What is even more important, the amplitude of the aperiodic oscillations 
grows rapidly. Therefore the finite-dimensional system is not applicable in this region and 
the chaotic solution shown in Fig.|6|b does not have a counterpart in the full nonreduced 
partial differential equation. 



7 Conclusions 



It was proposed by various authors [18, B3L 34 



specifically, to its long-wavelength component |33], [3 



, 41 1 that the coupling to the radiation (more 
5], j38| ) is the key factor determining 



the internal dynamics of the damped driven soliton, in particular its instability, bifurcations 
and transition to chaos. The main objective of the present work was to study the role of the 
soliton-radiation interaction within an approach which allows a rigorous decomposition of the 
phase space into the soliton and radiation modes. Our treatment is based on the Riemann- 
Hilbert problem, a modern version of the Inverse Scattering transform. Although the Inverse 
Scattering method had already been utilised in a related context (for the externally driven 
fill), our approach is different in that we are not assuming the damping and 



nls nil m 

driving terms to be small in any sense. Instead, we are exploiting a remarkable coincidence 
between the mathematical formula for the parametrically driven damped soliton and that 



of the soliton of the unperturbed NLS. This coincidence allows to associate the stationary 
damped-driven soliton with a stationary zero of the Riemann-Hilbert problem — for any h 
and 7. The evolution of all nearby solutions ip{ x ,t) can therefore be studied through the 
evolution of the corresponding scattering data. 

Some interesting insights are gained already from the linearised equations for the spectral 
data. Conclusions of this analysis are consistent with results of the linearisation in the 
space of fields ip(x,t) however the linearisation in the spectral space has an important 
advantage that it can be carried out analytically whereas linearised equations for ip(x,t) 
could only be studied by means of computer. The stationary soliton looses its stability as 
a result of its coupling to radiation waves. This had already been proposed before |I2| but 
now we have put this claim on rigorous footing. 

We attempted to advance beyond the linear approximation and track the effect of the 
long-wavelength radiation on the nonlinear dynamics of the soliton. To find a closed system 
for the evolution of the scattering data, we had to make two assumptions. First, we assumed 
that the solution of our damped driven NLS remains close to the stationary soliton of the 
same, perturbed, equation. Second, we assumed that the resulting system for the evolution 
of the spectral data is linear in radiation. (Note that we are not requiring the perturbation 
in the right-hand side of ([l]) to be small in any sense.) 

The outcome of this analysis is a dynamical system (|ST|)-(|&4]) comprising equations for the 
amplitude and phase of the soliton, and for the complex amplitude of the k = O-radiation. 
Finite-dimensional reductions of both externally [^, |34|, [35], 55] and parametrically H2|, [4~3 



driven NLS are available in literature and hence we need to emphasise the differences. The 
principal difference between our reduction and those obtained variationally [|4], [38], [42], [43| 



or by the Galerkin projections [35], lies in that our system of ODEs is not a product of 
any phenomenological Ansatz. It results from a rigorous expansion of the solution of the 
PDE over a set of "nonlinear modes" and then retaining only those modes whose effect 
on the dynamics we are trying to track down. A failure of one or the other variational or 
Galerkin approximation to capture essentials of the supercritical dynamics of the soliton 
does not provide any information on why this particular Ansatz fails. One has to try a 
variety of different Ansatze, select the one that gives the best fit and then attempt to make 
some semi- intuitive conclusions on the role of this or that ingredient of the trial function. 
On the contrary, our approach allows to explore, systematically, each part of the phase 
space and identify the nonlinear modes responsible for each particular dynamical effect. 



Next, our reduction technique is different from the approach of an influential paper |33 



which is also partially based on the Inverse Scattering. Besides the fact that the analysis 
of Ref. ]33 relies on the smallness of the perturbation, it only uses the Inverse Scattering 
transform to obtain the functional form of the radiation wave whereas its interaction with 
the soliton is introduced variationally. The damping term for the radiation was not part of 
the variational algorithm and had to be added in an ad-hoc way. The method does not define 
the amplitude of the radiation either; this is introduced empirically and then fitted to match 



the numerical data. Unlike Ref. [S3], our finite-dimensional reduction is uniquely defined by 
the choice of the ingredients of the localised attractor and does not require introduction of 
any phenomenological terms or fitting parameters. This uniqueness is reflected by the fact 
that (the linearisation of) our reduced system retains the self-similarity invariance of the 



(linearised) PDE. 

The analysis of the reduced dynamical system shows that it is capable of explaining only 
some parts and only some rough features of the attractor chart of the parametrically driven 



damped NLS []17 |. Most notably, the interaction with long radiation waves is sufficient to 
reproduce the approximate sequence of attractors arising when the driving strength is in- 
creased under the fixed dissipation coefficient. Consistently with computer simulations of 



17[ [39] , the finite-dimensional system exhibits the sequence "stationary soliton — > periodi- 
cally oscillating soliton" for larger dampings and "stationary soliton — > oscillating soliton — > 
flat zero solution" for smaller 7. For 7 = the reduced system does not predict oscillating 
solitons with bounded amplitudes (apart from a tiny window of h values). This is also in 
agreement with the behaviour observed in the full PDE [18]. However, the finite-dimensional 
system predicts — erroneously — the occurrence of the second, restabilising, Hopf bifurca- 
tion. As a result of that, the finite-dimensional fixed point turns out to be stable for all 7 
as long as h > 0.4, whereas the actual stability domain of the stationary soliton is much 
smaller One could have expected that the reduced and infinite-dimensional dynamics 
would be close in the region adjacent to the first, destabilising, Hopf bifurcation curve - 
which does provide a reasonably good approximation to the Hopf bifurcation curve in the 
full PDE. In the actual fact, however, details of the attractors and bifurcation sequences in 
the two systems are quite different even in that region. (An exception is the band of very 
small 7, 7 < 0.05.) 

Thus, taking into account just the k = component of radiation is insufficient for 
reproducing the entire complexity in the damped-driven soliton's dynamics. It is quite 
possible that the radiation waves with the frequency close to the double frequency of the 
soliton's linear oscillations [18| play a more important (or equally important) role than those 
with k = 0. On the other hand, numerical simulations on finite intervals reveal the excitation 



of radiations with several wavenumbers j|T[|. It is not unprobable that a similar wavenumber 
selection and competition occur on the infinite line. Finally, there are indications that the 
oscillating soliton of the externally driven NLS, is in fact a bound state of two solitons of 



the unperturbed, integrable, NLS equation [38|. A similar mechanism may operate in the 



parametrically driven case as well. Our approach allows to test all these possibilities; we are 
planning to do so in future publications. 
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9 Appendix: Discrete eigenvalues of the operator Q 



Here we demonstrate the existence of, and derive a lower bound for, discrete eigenvalues of 
the operator Q, equation (|6T|). The smallest eigenvalue, E±, can be sought for as a minimum 



of the corresponding Rayleigh quotient: 

£l=min I^. (97 ) 

*l>(k) J dkip 2 

In terms of (p — (1 + k 2 )ip, this can be rewritten as 

/ dkip | (i + k 2 )- 2 [(i + nk 2 ) 2 - h 2 ] + 2Hnf\ <p 

E\ = min T . — — — — „ — n , (98) 

V (k) J dk(l + k 2 Y 2 p 2 ' v ; 

where we have introduced a symmetric operator T: 

£<p(k) = (1 + k 2 y l B(l + k 2 )- l ip{k). (99) 



(Here £> is given by (|68|) .) 

For H = the operator £7 does not have discrete eigenvalues but as if grows, at least 
one discrete eigenvalue detaches from the continuum. Indeed, in the region \£\ < 1, \k\ < 1 
the kernel of the operator !F, 



sinh [f (k - £)] (1 + k 2 f/ 2 {l + £ 2 ) 3 / 2 : 

satisfies 



F(£,k) < A , , (Ifc- 1). (100) 
v ' y - 4smh7r v 7 v y 



For an arbitrary even function /i(fc) the inequality ( |100| ) yields 



dkd£ h{k)T{k, i)h(£) < ~ 4 I / < M> ( > I • ( ! 01 ) 




-l -l 



Choosing a suitable trial function ip = ipnik) in the quotient (^) and using the estimate 
(|101|) , one readily concludes that for any positive H the quotient can take values lying below 
1 — H 2 , the edge of the continuous spectrum. Hence the operator Q has at least one discrete 
eigenvalue for H > 0. 

Next, the norm of the operator fl9"9"|) is bounded: 

oo oo 

ll^ll 2 = f f dkd£F 2 (e,k) < -. 



Therefore T is a Fredholm operator and its eigenvalues £ m , T$ m = £ m 4>mi are bounded: 

^<ll^l| 2 <i (102) 



Returning to the operator Q (|67]), its eigenvalue E n (n = 1,2, ...) satisfies 



En 



(103) 



where ip n (k) is the associated eigenfunction: Q^ n = E n ip n . Using inequality ( |102[ ) in equation 
(|103|) we obtain a lower bound on E n : 



E n > 



J dkip n (l - H 2 + 2HQF) cp n 
Jd^I 



>(1-H 2 



HQ 
1-H 2 V 2 



(104) 



where (p n = (1 + fc 2 )^- 
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